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

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: }