ghc/rts/gmp/mpz/fib_ui.c

   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 }