Return to fib_ui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Initial revision

/* mpz_fib_ui(result, n) -- Set RESULT to the Nth Fibonacci number.

Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The GNU MP Library 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 Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

#include "gmp.h"
#include "gmp-impl.h"

/* This is fast, but could be made somewhat faster and neater.
   The timing is somewhat fluctuating for even/odd sizes because
   of the extra hair used to save variables and operations.  Here
   are a few things one might want to address:
     1. Avoid using 4 intermediate variables in mpz_fib_bigcase.
     2. Call mpn functions directly.  Straightforward for these functions.
     3. Merge the three functions into one.

Said by Kevin:
   Consider using the Lucas numbers L[n] as an auxiliary sequence, making
   it possible to do the "doubling" operation in mpz_fib_bigcase with two
   squares rather than two multiplies.  The formulas are a little more
   complicated, something like the following (untested).

       F[2n] = ((F[n]+L[n])^2 - 6*F[n]^2 - 4*(-1)^n) / 2
       L[2n] = 5*F[n]^2 + 2*(-1)^n

       F[2n+1] = (F[2n] + L[2n]) / 2
       L[2n+1] = (5*F[2n] + L[2n]) / 2

   The Lucas number that comes for free here could even be returned.

   Maybe there's formulas with two squares using just F[n], but I don't
   know of any.
*/

/* Determine the needed storage for Fib(n).  */
#define FIB_SIZE(n) (((mp_size_t) ((n)*0.695)) / BITS_PER_MP_LIMB + 2)

static void mpz_fib_bigcase _PROTO ((mpz_t, mpz_t, unsigned long int));
static void mpz_fib_basecase _PROTO ((mpz_t, mpz_t, unsigned long int));


#ifndef FIB_THRESHOLD
#define FIB_THRESHOLD 60
#endif

void
#if __STDC__
mpz_fib_ui (mpz_t r, unsigned long int n)
#else
mpz_fib_ui (r, n)
     mpz_t r;
     unsigned long int n;
#endif
{
  if (n == 0)
    mpz_set_ui (r, 0);
  else
    {
      mpz_t t1;
      mpz_init (t1);
      if (n < FIB_THRESHOLD)
	mpz_fib_basecase (t1, r, n);
      else
	mpz_fib_bigcase (t1, r, n);
      mpz_clear (t1);
    }
}

static void
#if __STDC__
mpz_fib_basecase (mpz_t t1, mpz_t t2, unsigned long int n)
#else
mpz_fib_basecase (t1, t2, n)
     mpz_t t1;
     mpz_t t2;
     unsigned long int n;
#endif
{
  unsigned long int m, i;

  mpz_set_ui (t1, 0);
  mpz_set_ui (t2, 1);
  m = n/2;
  for (i = 0; i < m; i++)
    {
      mpz_add (t1, t1, t2);
      mpz_add (t2, t1, t2);
    }
  if ((n & 1) == 0)
    {
      mpz_sub (t1, t2, t1);
      mpz_sub (t2, t2, t1);	/* trick: recover t1 value just overwritten */
    }
}

static void
#if __STDC__
mpz_fib_bigcase (mpz_t t1, mpz_t t2, unsigned long int n)
#else
mpz_fib_bigcase (t1, t2, n)
     mpz_t t1;
     mpz_t t2;
     unsigned long int n;
#endif
{
  unsigned long int n2;
  int ni, i;
  mpz_t x1, x2, u1, u2;

  ni = 0;
  for (n2 = n; n2 >= FIB_THRESHOLD; n2 /= 2)
    ni++;

  mpz_fib_basecase (t1, t2, n2);

  mpz_init (x1);
  mpz_init (x2);
  mpz_init (u1);
  mpz_init (u2);

  for (i = ni - 1; i >= 0; i--)
    {
      mpz_mul_2exp (x1, t1, 1);
      mpz_mul_2exp (x2, t2, 1);

      mpz_add (x1, x1, t2);
      mpz_sub (x2, x2, t1);

      mpz_mul (u1, t2, x1);
      mpz_mul (u2, t1, x2);

      if (((n >> i) & 1) == 0)
	{
	  mpz_sub (t1, u1, u2);
	  mpz_set (t2, u1);
	}
      else
	{
	  mpz_set (t1, u1);
	  mpz_mul_2exp (t2, u1, 1);
	  mpz_sub (t2, t2, u2);
	}
    }

  mpz_clear (x1);
  mpz_clear (x2);
  mpz_clear (u1);
  mpz_clear (u2);
}