From ab9168949ec2762698d6293adf17b637f30b891e Mon Sep 17 00:00:00 2001 From: Mark Wooding Date: Sun, 7 Apr 2013 23:44:34 +0100 Subject: [PATCH] New function and example program computes Fibonacci numbers fairly fast. --- Makefile.m4 | 8 +- fibonacci.c | 145 ++++++++++++++++++++++++++++++ mp-fibonacci.c | 271 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ mp-fibonacci.h | 58 ++++++++++++ tests/mp | 18 ++++ 5 files changed, 497 insertions(+), 3 deletions(-) create mode 100644 fibonacci.c create mode 100644 mp-fibonacci.c create mode 100644 mp-fibonacci.h diff --git a/Makefile.m4 b/Makefile.m4 index c191886..d05e6ee 100644 --- a/Makefile.m4 +++ b/Makefile.m4 @@ -173,7 +173,7 @@ pkginclude_HEADERS = \ key.h key-error.h key-data.h passphrase.h pixie.h lmem.h \ mpx.h bitops.h mpw.h mpscan.h mparena.h mp.h mptext.h mpint.h \ exp.h mpbarrett.h mpmont.h mpreduce.h mp-exp.h \ - mpcrt.h mprand.h mpmul.h \ + mpcrt.h mprand.h mpmul.h mp-fibonacci.h \ gfx.h gf.h gfreduce.h gfn.h gf-exp.h \ primetab.h wheel.h pfilt.h rabin.h \ pgen.h primeiter.h prim.h strongprime.h limlee.h keycheck.h \ @@ -200,7 +200,7 @@ define(`MP_SOURCES', mp-sqrt.c mp-gcd.c mp-jacobi.c mp-modsqrt.c mp-exp.c mp-modexp.c \ mpint.c mptext-file.c mptext-dstr.c \ mptext-len.c \ - exp.c mpcrt.c mpmul.c mprand.c \ + exp.c mpcrt.c mpmul.c mprand.c mp-fibonacci.c \ mpbarrett.c mpbarrett-exp.c mpbarrett-mexp.c mpbarrett-exp.h \ mpmont.c mpmont-exp.c mpmont-mexp.c mpmont-exp.h \ mpreduce.c mpreduce-exp.h \ @@ -273,7 +273,7 @@ gfx-sqr.lo: gfx-sqr-tab.h ## --- Utility programs --- bin_PROGRAMS = \ - dsig key pixie cookie rspit factorial hashsum mkphrase \ + dsig key pixie cookie rspit factorial fibonacci hashsum mkphrase \ catcrypt catsign noinst_LIBRARIES = libcatcrypt.a pkgconfigdir = $(libdir)/pkgconfig @@ -302,6 +302,7 @@ key_SOURCES = keyutil.c hashsum_SOURCES = hashsum.c rspit_SOURCES = rspit.c factorial_SOURCES = factorial.c +fibonacci_SOURCES = fibonacci.c perftest_SOURCES = perftest.c perftest_LDADD = $(CATACOMB_LIBS) $(LDADD) pixie_SOURCES = pixie.c pixie-common.c lmem.c arena.c passphrase.c @@ -440,6 +441,7 @@ CTESTRIG(mpmont-mexp) CTESTRIG(mpreduce) CTESTRIG(mpcrt) CTESTRIG(mpmul) +CTESTRIG(mp-fibonacci) CTESTRIG(rsa-test) CTESTRIG(gfx) CTESTRIG(gfx-sqr) diff --git a/fibonacci.c b/fibonacci.c new file mode 100644 index 0000000..4b5a958 --- /dev/null +++ b/fibonacci.c @@ -0,0 +1,145 @@ +/* -*-c-*- + * + * Example Fibonacci number computation + * + * (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 "config.h" + +#include +#include +#include +#include + +#include +#include +#include + +#include "mpint.h" +#include "mptext.h" + +#include "mp-fibonacci.h" + +/*----- Main code ---------------------------------------------------------*/ + +static void usage(FILE *fp) + { pquis(fp, "Usage: $ [-r RADIX] INTEGER\n"); } + +static void version(FILE *fp) + { pquis(fp, "$, Catacomb version " VERSION "\n"); } + +static void help(FILE *fp) +{ + version(fp); + putc('\n', fp); + usage(fp); + fputs("\n\ +Prints the Nth Fibonacci number for a given inteer N. Input may be in\n\ +decimal (the default), octal with preceding zero, hex with preceding `0x',\n\ +or any base N between 2 and 62 inclusive with preceding `N_'. Output\n\ +may be in any base between 2 and 62; the default is base 10. For bases\n\ +between 11 and 36 inclusive, lowercase letters of either case are used as\n\ +additional digits with values 10 upwards; lowercase is always used for\n\ +output. For bases between 37 and 62 inclusive, lowercase letters have\n\ +lower value than uppercase; hence `a' has the value 10, while `A' has the\n\ +value 36.\n\ +\n\ +Options provided:\n\ +\n\ +-h, --help Display this help message.\n\ +-v, --version Display the version number.\n\ +-u, --usage Display a usage message.\n\ +\n\ +-r, --radix=N Write output in base N.\n\ +", fp); +} + +int main(int argc, char *argv[]) +{ + long n; + int r = 10; + char *p; + mp *f, *lmin, *lmax, *nn; + unsigned fl = 0; + +#define f_bogus 1u + + ego(argv[0]); + + for (;;) { + static const struct option opt[] = { + { "help", 0, 0, 'h' }, + { "version", 0, 0, 'v' }, + { "usage", 0, 0, 'u' }, + { "radix", OPTF_ARGREQ, 0, 'r' }, + { 0, 0, 0, 0 } + }; + int i = mdwopt(argc, argv, "hvur:", opt, 0, 0, 0); + if (i < 0) + break; + switch (i) { + case 'h': + help(stdout); + exit(0); + case 'v': + version(stdout); + exit(0); + case 'u': + usage(stdout); + exit(0); + case 'r': + r = atoi(optarg); + if (r < 2 || r > 62) + die(EXIT_FAILURE, "bad radix `%s'", optarg); + break; + default: + fl |= f_bogus; + break; + } + } + + if (optind + 1 != argc || (fl & f_bogus)) { + usage(stderr); + exit(EXIT_FAILURE); + } + lmin = mp_fromlong(MP_NEW, LONG_MIN); + lmax = mp_fromlong(MP_NEW, LONG_MAX); + p = argv[optind]; + while (isspace((unsigned char)*p)) p++; + nn = mp_readstring(MP_NEW, argv[optind], &p, 0); + while (isspace((unsigned char)*p)) p++; + if (!nn || *p || MP_CMP(lmin, >, nn) || MP_CMP(nn, >, lmax)) + die(EXIT_FAILURE, "bad integer `%s'", argv[optind]); + n = mp_tolong(nn); + mp_drop(nn); mp_drop(lmin); mp_drop(lmax); + f = mp_fibonacci(n); + mp_writefile(f, stdout, r); + fputc('\n', stdout); + mp_drop(f); + return (0); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-fibonacci.c b/mp-fibonacci.c new file mode 100644 index 0000000..e103979 --- /dev/null +++ b/mp-fibonacci.c @@ -0,0 +1,271 @@ +/* -*-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 + +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 -------------------------------------------------*/ diff --git a/mp-fibonacci.h b/mp-fibonacci.h new file mode 100644 index 0000000..bdc20a6 --- /dev/null +++ b/mp-fibonacci.h @@ -0,0 +1,58 @@ +/* -*-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. + */ + +#ifndef CATACOMB_MP_FIBONACCI_H +#define CATACOMB_MP_FIBONACCI_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef CATACOMB_MP_H +# include "mp.h" +#endif + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @mp_fibonacci@ --- * + * + * Arguments: @long n@ = index desired (may be negative) + * + * Returns: The %$n$%th Fibonacci number. + */ + +extern mp *mp_fibonacci(long /*n*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/tests/mp b/tests/mp index 1800fea..0fbaf06 100644 --- a/tests/mp +++ b/tests/mp @@ -308,3 +308,21 @@ factorial { 500 1220136825991110068701238785423046926253574342803192842192413588385845373153881997605496447502203281863013616477148203584163378722078177200480785205159329285477907571939330603772960859086270429174547882424912726344305670173270769461062802310452644218878789465754777149863494367781037644274033827365397471386477878495438489595537537990423241061271326984327745715546309977202781014561081188373709531016356324432987029563896628911658974769572087926928871281780070265174507768410719624390394322536422605234945850129918571501248706961568141625359056693423813008856249246891564126775654481886506593847951775360894005745238940335798476363944905313062323749066445048824665075946735862074637925184200459369692981022263971952597190945217823331756934581508552332820762820023402626907898342451712006207714640979456116127629145951237229913340169552363850942885592018727433795173014586357570828355780158735432768888680120399882384702151467605445407663535984174430480128938313896881639487469658817504506926365338175055478128640000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000; } + +fibonacci { + -20 -6765; + -19 4181; + -10 -55; + -9 34; + -2 -1; + -1 1; + 0 0; + 1 1; + 2 1; + 5 5; + 10 55; + 19 4181; + 20 6765; + 100 354224848179261915075; + 1000 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875; +} -- 2.11.0