contrib/gmp/mpn/generic/fib2_ui.c

   1 /* mpn_fib2_ui -- calculate Fibonacci numbers.
   2
   3    THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
   4    CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
   5    FUTURE GNU MP RELEASES.
   6
   7 Copyright 2001, 2002, 2005 Free Software Foundation, Inc.
   8
   9 This file is part of the GNU MP Library.
  10
  11 The GNU MP Library is free software; you can redistribute it and/or modify
  12 it under the terms of the GNU Lesser General Public License as published by
  13 the Free Software Foundation; either version 3 of the License, or (at your
  14 option) any later version.
  15
  16 The GNU MP Library is distributed in the hope that it will be useful, but
  17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  19 License for more details.
  20
  21 You should have received a copy of the GNU Lesser General Public License
  22 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  23
  24 #include <stdio.h>
  25 #include "gmp.h"
  26 #include "gmp-impl.h"
  27 #include "longlong.h"
  28
  29
  30 /* change this to "#define TRACE(x) x" for diagnostics */
  31 #define TRACE(x)
  32
  33
  34 /* Store F[n] at fp and F[n-1] at f1p.  fp and f1p should have room for
  35    MPN_FIB2_SIZE(n) limbs.
  36
  37    The return value is the actual number of limbs stored, this will be at
  38    least 1.  fp[size-1] will be non-zero, except when n==0, in which case
  39    fp[0] is 0 and f1p[0] is 1.  f1p[size-1] can be zero, since F[n-1]<F[n]
  40    (for n>0).
  41
  42    Notes:
  43
  44    In F[2k+1] with k even, +2 is applied to 4*F[k]^2 just by ORing into the
  45    low limb.
  46
  47    In F[2k+1] with k odd, -2 is applied to the low limb of 4*F[k]^2 -
  48    F[k-1]^2.  This F[2k+1] is an F[4m+3] and such numbers are congruent to
  49    1, 2 or 5 mod 8, which means no underflow reaching it with a -2 (since
  50    that would leave 6 or 7 mod 8).
  51
  52    This property of F[4m+3] can be verified by induction on F[4m+3] =
  53    7*F[4m-1] - F[4m-5], that formula being a standard lucas sequence
  54    identity U[i+j] = U[i]*V[j] - U[i-j]*Q^j.
  55
  56    Enhancements:
  57
  58    If there was an mpn_addlshift, it'd be possible to eliminate the yp
  59    temporary, using xp=F[k]^2, fp=F[k-1]^2, f1p=xp+fp, fp+=4*fp, fp-=f1p,
  60    fp+=2*(-1)^n, etc.  */
  61
  62 mp_size_t
  63 mpn_fib2_ui (mp_ptr fp, mp_ptr f1p, unsigned long int n)
  64 {
  65   mp_ptr         xp, yp;
  66   mp_size_t      size;
  67   unsigned long  nfirst, mask;
  68   TMP_DECL;
  69
  70   TRACE (printf ("mpn_fib2_ui n=%lu\n", n));
  71
  72   ASSERT (! MPN_OVERLAP_P (fp, MPN_FIB2_SIZE(n), f1p, MPN_FIB2_SIZE(n)));
  73
  74   /* Take a starting pair from the table. */
  75   mask = 1;
  76   for (nfirst = n; nfirst > FIB_TABLE_LIMIT; nfirst /= 2)
  77     mask <<= 1;
  78   TRACE (printf ("nfirst=%lu mask=0x%lX\n", nfirst, mask));
  79
  80   f1p[0] = FIB_TABLE ((int) nfirst - 1);
  81   fp[0]  = FIB_TABLE (nfirst);
  82   size = 1;
  83
  84   /* Skip to the end if the table lookup gives the final answer. */
  85   if (mask != 1)
  86     {
  87       mp_size_t  alloc;
  88
  89       TMP_MARK;
  90       alloc = MPN_FIB2_SIZE (n);
  91       TMP_ALLOC_LIMBS_2 (xp,alloc, yp,alloc);
  92
  93       do
  94         {
  95           mp_limb_t  c;
  96
  97           /* Here fp==F[k] and f1p==F[k-1], with k being the bits of n from
  98              n&mask upwards.
  99
 100              The next bit of n is n&(mask>>1) and we'll double to the pair
 101              fp==F[2k],f1p==F[2k-1] or fp==F[2k+1],f1p==F[2k], according as
 102              that bit is 0 or 1 respectively.  */
 103
 104           TRACE (printf ("k=%lu mask=0x%lX size=%ld alloc=%ld\n",
 105                          n >> refmpn_count_trailing_zeros(mask),
 106                          mask, size, alloc);
 107                  mpn_trace ("fp ", fp, size);
 108                  mpn_trace ("f1p", f1p, size));
 109
 110           /* fp normalized, f1p at most one high zero */
 111           ASSERT (fp[size-1] != 0);
 112           ASSERT (f1p[size-1] != 0 || f1p[size-2] != 0);
 113
 114           /* f1p[size-1] might be zero, but this occurs rarely, so it's not
 115              worth bothering checking for it */
 116           ASSERT (alloc >= 2*size);
 117           mpn_sqr_n (xp, fp,  size);
 118           mpn_sqr_n (yp, f1p, size);
 119           size *= 2;
 120
 121           /* Shrink if possible.  Since fp was normalized there'll be at
 122              most one high zero on xp (and if there is then there's one on
 123              yp too).  */
 124           ASSERT (xp[size-1] != 0 || yp[size-1] == 0);
 125           size -= (xp[size-1] == 0);
 126           ASSERT (xp[size-1] != 0);  /* only one xp high zero */
 127
 128           /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k.
 129              n&mask is the low bit of our implied k.  */
 130           c = mpn_lshift (fp, xp, size, 2);
 131           fp[0] |= (n & mask ? 0 : 2);   /* possible +2 */
 132           c -= mpn_sub_n (fp, fp, yp, size);
 133           ASSERT (n & (mask << 1) ? fp[0] != 0 && fp[0] != 1 : 1);
 134           fp[0] -= (n & mask ? 2 : 0);   /* possible -2 */
 135           ASSERT (alloc >= size+1);
 136           xp[size] = 0;
 137           yp[size] = 0;
 138           fp[size] = c;
 139           size += (c != 0);
 140
 141           /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2.
 142              F[2k-1]<F[2k+1] so no carry out of "size" limbs. */
 143           ASSERT_NOCARRY (mpn_add_n (f1p, xp, yp, size));
 144
 145           /* now n&mask is the new bit of n being considered */
 146           mask >>= 1;
 147
 148           /* Calculate F[2k] = F[2k+1] - F[2k-1], replacing the unwanted one of
 149              F[2k+1] and F[2k-1].  */
 150           ASSERT_NOCARRY (mpn_sub_n ((n & mask ? f1p : fp), fp, f1p, size));
 151
 152           /* Can have a high zero after replacing F[2k+1] with F[2k].
 153              f1p will have a high zero if fp does. */
 154           ASSERT (fp[size-1] != 0 || f1p[size-1] == 0);
 155           size -= (fp[size-1] == 0);
 156         }
 157       while (mask != 1);
 158
 159       TMP_FREE;
 160     }
 161
 162   TRACE (printf ("done size=%ld\n", size);
 163          mpn_trace ("fp ", fp, size);
 164          mpn_trace ("f1p", f1p, size));
 165
 166   return size;
 167 }