Return to fib_ui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/fib_ui.c, Revision 1.1

1.1     ! maekawa     1: /* mpz_fib_ui(result, n) -- Set RESULT to the Nth Fibonacci number.
        !             2:
        !             3: Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.
        !             4:
        !             5: This file is part of the GNU MP Library.
        !             6:
        !             7: The GNU MP Library is free software; you can redistribute it and/or modify
        !             8: it under the terms of the GNU Lesser General Public License as published by
        !             9: the Free Software Foundation; either version 2.1 of the License, or (at your
        !            10: option) any later version.
        !            11:
        !            12: The GNU MP Library is distributed in the hope that it will be useful, but
        !            13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
        !            14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
        !            15: License for more details.
        !            16:
        !            17: You should have received a copy of the GNU Lesser General Public License
        !            18: along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
        !            19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
        !            20: MA 02111-1307, USA. */
        !            21:
        !            22: #include "gmp.h"
        !            23: #include "gmp-impl.h"
        !            24:
        !            25: /* This is fast, but could be made somewhat faster and neater.
        !            26:    The timing is somewhat fluctuating for even/odd sizes because
        !            27:    of the extra hair used to save variables and operations.  Here
        !            28:    are a few things one might want to address:
        !            29:      1. Avoid using 4 intermediate variables in mpz_fib_bigcase.
        !            30:      2. Call mpn functions directly.  Straightforward for these functions.
        !            31:      3. Merge the three functions into one.
        !            32:
        !            33: Said by Kevin:
        !            34:    Consider using the Lucas numbers L[n] as an auxiliary sequence, making
        !            35:    it possible to do the "doubling" operation in mpz_fib_bigcase with two
        !            36:    squares rather than two multiplies.  The formulas are a little more
        !            37:    complicated, something like the following (untested).
        !            38:
        !            39:        F[2n] = ((F[n]+L[n])^2 - 6*F[n]^2 - 4*(-1)^n) / 2
        !            40:        L[2n] = 5*F[n]^2 + 2*(-1)^n
        !            41:
        !            42:        F[2n+1] = (F[2n] + L[2n]) / 2
        !            43:        L[2n+1] = (5*F[2n] + L[2n]) / 2
        !            44:
        !            45:    The Lucas number that comes for free here could even be returned.
        !            46:
        !            47:    Maybe there's formulas with two squares using just F[n], but I don't
        !            48:    know of any.
        !            49: */
        !            50:
        !            51: /* Determine the needed storage for Fib(n).  */
        !            52: #define FIB_SIZE(n) (((mp_size_t) ((n)*0.695)) / BITS_PER_MP_LIMB + 2)
        !            53:
        !            54: static void mpz_fib_bigcase _PROTO ((mpz_t, mpz_t, unsigned long int));
        !            55: static void mpz_fib_basecase _PROTO ((mpz_t, mpz_t, unsigned long int));
        !            56:
        !            57:
        !            58: #ifndef FIB_THRESHOLD
        !            59: #define FIB_THRESHOLD 60
        !            60: #endif
        !            61:
        !            62: void
        !            63: #if __STDC__
        !            64: mpz_fib_ui (mpz_t r, unsigned long int n)
        !            65: #else
        !            66: mpz_fib_ui (r, n)
        !            67:      mpz_t r;
        !            68:      unsigned long int n;
        !            69: #endif
        !            70: {
        !            71:   if (n == 0)
        !            72:     mpz_set_ui (r, 0);
        !            73:   else
        !            74:     {
        !            75:       mpz_t t1;
        !            76:       mpz_init (t1);
        !            77:       if (n < FIB_THRESHOLD)
        !            78:        mpz_fib_basecase (t1, r, n);
        !            79:       else
        !            80:        mpz_fib_bigcase (t1, r, n);
        !            81:       mpz_clear (t1);
        !            82:     }
        !            83: }
        !            84:
        !            85: static void
        !            86: #if __STDC__
        !            87: mpz_fib_basecase (mpz_t t1, mpz_t t2, unsigned long int n)
        !            88: #else
        !            89: mpz_fib_basecase (t1, t2, n)
        !            90:      mpz_t t1;
        !            91:      mpz_t t2;
        !            92:      unsigned long int n;
        !            93: #endif
        !            94: {
        !            95:   unsigned long int m, i;
        !            96:
        !            97:   mpz_set_ui (t1, 0);
        !            98:   mpz_set_ui (t2, 1);
        !            99:   m = n/2;
        !           100:   for (i = 0; i < m; i++)
        !           101:     {
        !           102:       mpz_add (t1, t1, t2);
        !           103:       mpz_add (t2, t1, t2);
        !           104:     }
        !           105:   if ((n & 1) == 0)
        !           106:     {
        !           107:       mpz_sub (t1, t2, t1);
        !           108:       mpz_sub (t2, t2, t1);    /* trick: recover t1 value just overwritten */
        !           109:     }
        !           110: }
        !           111:
        !           112: static void
        !           113: #if __STDC__
        !           114: mpz_fib_bigcase (mpz_t t1, mpz_t t2, unsigned long int n)
        !           115: #else
        !           116: mpz_fib_bigcase (t1, t2, n)
        !           117:      mpz_t t1;
        !           118:      mpz_t t2;
        !           119:      unsigned long int n;
        !           120: #endif
        !           121: {
        !           122:   unsigned long int n2;
        !           123:   int ni, i;
        !           124:   mpz_t x1, x2, u1, u2;
        !           125:
        !           126:   ni = 0;
        !           127:   for (n2 = n; n2 >= FIB_THRESHOLD; n2 /= 2)
        !           128:     ni++;
        !           129:
        !           130:   mpz_fib_basecase (t1, t2, n2);
        !           131:
        !           132:   mpz_init (x1);
        !           133:   mpz_init (x2);
        !           134:   mpz_init (u1);
        !           135:   mpz_init (u2);
        !           136:
        !           137:   for (i = ni - 1; i >= 0; i--)
        !           138:     {
        !           139:       mpz_mul_2exp (x1, t1, 1);
        !           140:       mpz_mul_2exp (x2, t2, 1);
        !           141:
        !           142:       mpz_add (x1, x1, t2);
        !           143:       mpz_sub (x2, x2, t1);
        !           144:
        !           145:       mpz_mul (u1, t2, x1);
        !           146:       mpz_mul (u2, t1, x2);
        !           147:
        !           148:       if (((n >> i) & 1) == 0)
        !           149:        {
        !           150:          mpz_sub (t1, u1, u2);
        !           151:          mpz_set (t2, u1);
        !           152:        }
        !           153:       else
        !           154:        {
        !           155:          mpz_set (t1, u1);
        !           156:          mpz_mul_2exp (t2, u1, 1);
        !           157:          mpz_sub (t2, t2, u2);
        !           158:        }
        !           159:     }
        !           160:
        !           161:   mpz_clear (x1);
        !           162:   mpz_clear (x2);
        !           163:   mpz_clear (u1);
        !           164:   mpz_clear (u2);
        !           165: }