fib_ui.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* mpz_fib_ui -- calculate Fibonacci numbers.
  2. Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc.
  3. This file is part of the GNU MP Library.
  4. The GNU MP Library is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU Lesser General Public License as published by
  6. the Free Software Foundation; either version 3 of the License, or (at your
  7. option) any later version.
  8. The GNU MP Library is distributed in the hope that it will be useful, but
  9. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. License for more details.
  12. You should have received a copy of the GNU Lesser General Public License
  13. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  14. #include <stdio.h>
  15. #include "gmp.h"
  16. #include "gmp-impl.h"
  17. #include "longlong.h"
  18. /* change to "#define TRACE(x) x" to get some traces */
  19. #define TRACE(x)
  20. /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low
  21.    limb because the result F[2k+1] is an F[4m+3] and such numbers are always
  22.    == 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7.  (This is
  23.    the same as in mpn_fib2_ui.)
  24.    In the F[2k+1] for k even, the +2 won't give a carry out of the low limb
  25.    in normal circumstances.  This is an F[4m+1] and we claim that F[3*2^b+1]
  26.    == 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence
  27.    if n < 2^GMP_NUMB_BITS then F[n] cannot have a low limb of 0 or 1.  No
  28.    proof for this claim, but it's been verified up to b==32 and has such a
  29.    nice pattern it must be true :-).  Of interest is that F[3*2^b] == 0 mod
  30.    2^(b+1) seems to hold too.
  31.    When n >= 2^GMP_NUMB_BITS, which can arise in a nails build, then the low
  32.    limb of F[4m+1] can certainly be 1, and an mpn_add_1 must be used.  */
  33. void
  34. mpz_fib_ui (mpz_ptr fn, unsigned long n)
  35. {
  36.   mp_ptr         fp, xp, yp;
  37.   mp_size_t      size, xalloc;
  38.   unsigned long  n2;
  39.   mp_limb_t      c, c2;
  40.   TMP_DECL;
  41.   if (n <= FIB_TABLE_LIMIT)
  42.     {
  43.       PTR(fn)[0] = FIB_TABLE (n);
  44.       SIZ(fn) = (n != 0);      /* F[0]==0, others are !=0 */
  45.       return;
  46.     }
  47.   n2 = n/2;
  48.   xalloc = MPN_FIB2_SIZE (n2) + 1;
  49.   MPZ_REALLOC (fn, 2*xalloc+1);
  50.   fp = PTR (fn);
  51.   TMP_MARK;
  52.   TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc);
  53.   size = mpn_fib2_ui (xp, yp, n2);
  54.   TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lun",
  55.                  n >> 1, size, n&1);
  56.          mpn_trace ("xp", xp, size);
  57.          mpn_trace ("yp", yp, size));
  58.   if (n & 1)
  59.     {
  60.       /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k  */
  61.       mp_size_t  xsize, ysize;
  62. #if HAVE_NATIVE_mpn_add_n_sub_n
  63.       xp[size] = mpn_lshift (xp, xp, size, 1);
  64.       yp[size] = 0;
  65.       ASSERT_NOCARRY (mpn_add_n_sub_n (xp, yp, xp, yp, size+1));
  66.       xsize = size + (xp[size] != 0);
  67.       ysize = size + (yp[size] != 0);
  68. #else
  69.       c2 = mpn_lshift (fp, xp, size, 1);
  70.       c = c2 + mpn_add_n (xp, fp, yp, size);
  71.       xp[size] = c;
  72.       xsize = size + (c != 0);
  73.       c2 -= mpn_sub_n (yp, fp, yp, size);
  74.       yp[size] = c2;
  75.       ASSERT (c2 <= 1);
  76.       ysize = size + c2;
  77. #endif
  78.       size = xsize + ysize;
  79.       c = mpn_mul (fp, xp, xsize, yp, ysize);
  80. #if GMP_NUMB_BITS >= BITS_PER_ULONG
  81.       /* no overflow, see comments above */
  82.       ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= GMP_NUMB_MAX-2);
  83.       fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2));
  84. #else
  85.       if (n & 2)
  86.         {
  87.           ASSERT (fp[0] >= 2);
  88.           fp[0] -= 2;
  89.         }
  90.       else
  91.         {
  92.           ASSERT (c != GMP_NUMB_MAX); /* because it's the high of a mul */
  93.           c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2));
  94.           fp[size-1] = c;
  95.         }
  96. #endif
  97.     }
  98.   else
  99.     {
  100.       /* F[2k] = F[k]*(F[k]+2F[k-1]) */
  101.       mp_size_t  xsize, ysize;
  102.       c = mpn_lshift (yp, yp, size, 1);
  103.       c += mpn_add_n (yp, yp, xp, size);
  104.       yp[size] = c;
  105.       xsize = size;
  106.       ysize = size + (c != 0);
  107.       size += ysize;
  108.       c = mpn_mul (fp, yp, ysize, xp, xsize);
  109.     }
  110.   /* one or two high zeros */
  111.   size -= (c == 0);
  112.   size -= (fp[size-1] == 0);
  113.   SIZ(fn) = size;
  114.   TRACE (printf ("done special, size=%ldn", size);
  115.          mpn_trace ("fp ", fp, size));
  116.   TMP_FREE;
  117. }