Annotation of OpenXM_contrib/gmp/mpz/fib_ui.c, Revision 1.1.1.2
1.1.1.2 ! ohara 1: /* mpz_fib_ui -- calculate Fibonacci numbers.
1.1 maekawa 2:
1.1.1.2 ! ohara 3: Copyright 2000, 2001 Free Software Foundation, Inc.
1.1 maekawa 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:
1.1.1.2 ! ohara 22: #include <stdio.h>
1.1 maekawa 23: #include "gmp.h"
24: #include "gmp-impl.h"
1.1.1.2 ! ohara 25: #include "longlong.h"
1.1 maekawa 26:
27:
1.1.1.2 ! ohara 28: /* change to "#define TRACE(x) x" to get some traces */
! 29: #define TRACE(x)
1.1 maekawa 30:
31:
1.1.1.2 ! ohara 32: /* In the F[2k+1] below for k odd, the -2 won't give a borrow from the low
! 33: limb because the result F[2k+1] is an F[4m+3] and such numbers are always
! 34: == 1, 2 or 5 mod 8, whereas an underflow would leave 6 or 7. (This is
! 35: the same as in mpn_fib2_ui.)
1.1 maekawa 36:
1.1.1.2 ! ohara 37: In the F[2k+1] for k even, the +2 won't give a carry out of the low limb
! 38: in normal circumstances. This is an F[4m+1] and we claim that F[3*2^b+1]
! 39: == 1 mod 2^b is the first F[4m+1] congruent to 0 or 1 mod 2^b, and hence
! 40: if n < 2^BITS_PER_MP_LIMB then F[n] cannot have a low limb of 0 or 1. No
! 41: proof for this claim, but it's been verified up to b==32 and has such a
! 42: nice pattern it must be true :-). Of interest is that F[3*2^b] == 0 mod
! 43: 2^(b+1) seems to hold too.
1.1 maekawa 44:
1.1.1.2 ! ohara 45: When n >= 2^BITS_PER_MP_LIMB, which can arise in test setups with a small
! 46: limb, then the low limb of F[4m+1] can certainly be 1, and an mpn_add_1
! 47: must be used. */
1.1 maekawa 48:
49: void
1.1.1.2 ! ohara 50: mpz_fib_ui (mpz_ptr fn, unsigned long n)
1.1 maekawa 51: {
1.1.1.2 ! ohara 52: mp_ptr fp, xp, yp;
! 53: mp_size_t size, xalloc;
! 54: unsigned long n2;
! 55: mp_limb_t c, c2;
! 56: TMP_DECL (marker);
! 57:
! 58: if (n <= FIB_TABLE_LIMIT)
1.1 maekawa 59: {
1.1.1.2 ! ohara 60: PTR(fn)[0] = FIB_TABLE (n);
! 61: SIZ(fn) = (n != 0); /* F[0]==0, others are !=0 */
! 62: return;
1.1 maekawa 63: }
64:
1.1.1.2 ! ohara 65: n2 = n/2;
! 66: xalloc = MPN_FIB2_SIZE (n2) + 1;
! 67: MPZ_REALLOC (fn, 2*xalloc+1);
! 68: fp = PTR (fn);
! 69:
! 70: TMP_MARK (marker);
! 71: TMP_ALLOC_LIMBS_2 (xp,xalloc, yp,xalloc);
! 72: size = mpn_fib2_ui (xp, yp, n2);
! 73:
! 74: TRACE (printf ("mpz_fib_ui last step n=%lu size=%ld bit=%lu\n",
! 75: n >> 1, size, n&1);
! 76: mpn_trace ("xp", xp, size);
! 77: mpn_trace ("yp", yp, size));
1.1 maekawa 78:
1.1.1.2 ! ohara 79: if (n & 1)
1.1 maekawa 80: {
1.1.1.2 ! ohara 81: /* F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k */
! 82: mp_size_t xsize, ysize;
1.1 maekawa 83:
1.1.1.2 ! ohara 84: #if HAVE_NATIVE_mpn_addsub_n
! 85: xp[size] = mpn_lshift (xp, xp, size, 1);
! 86: yp[size] = 0;
! 87: ASSERT_NOCARRY (mpn_addsub_n (xp, yp, xp, yp, size+1));
! 88: xsize = size + (xp[size] != 0);
! 89: ysize = size + (yp[size] != 0);
1.1 maekawa 90: #else
1.1.1.2 ! ohara 91: c2 = mpn_lshift (fp, xp, size, 1);
! 92: c = c2 + mpn_add_n (xp, fp, yp, size);
! 93: xp[size] = c;
! 94: xsize = size + (c != 0);
! 95: c2 -= mpn_sub_n (yp, fp, yp, size);
! 96: yp[size] = c2;
! 97: ASSERT (c2 <= 1);
! 98: ysize = size + c2;
1.1 maekawa 99: #endif
100:
1.1.1.2 ! ohara 101: size = xsize + ysize;
! 102: c = mpn_mul (fp, xp, xsize, yp, ysize);
! 103:
! 104: #if BITS_PER_MP_LIMB >= BITS_PER_ULONG
! 105: /* no overflow, see comments above */
! 106: ASSERT (n & 2 ? fp[0] >= 2 : fp[0] <= MP_LIMB_T_MAX-2);
! 107: fp[0] += (n & 2 ? -CNST_LIMB(2) : CNST_LIMB(2));
! 108: #else
! 109: /* this code only for testing with small limbs, limb<ulong is unusual */
! 110: if (n & 2)
! 111: {
! 112: ASSERT (fp[0] >= 2);
! 113: fp[0] -= 2;
! 114: }
! 115: else
! 116: {
! 117: ASSERT (c != MP_LIMB_T_MAX); /* because it's the high of a mul */
! 118: c += mpn_add_1 (fp, fp, size-1, CNST_LIMB(2));
! 119: fp[size-1] = c;
! 120: }
! 121: #endif
! 122: }
! 123: else
1.1 maekawa 124: {
1.1.1.2 ! ohara 125: /* F[2k] = F[k]*(F[k]+2F[k-1]) */
1.1 maekawa 126:
1.1.1.2 ! ohara 127: mp_size_t xsize, ysize;
! 128: c = mpn_lshift (yp, yp, size, 1);
! 129: c += mpn_add_n (yp, yp, xp, size);
! 130: yp[size] = c;
! 131: xsize = size;
! 132: ysize = size + (c != 0);
! 133: size += ysize;
! 134: c = mpn_mul (fp, yp, ysize, xp, xsize);
! 135: }
1.1 maekawa 136:
1.1.1.2 ! ohara 137: /* one or two high zeros */
! 138: size -= (c == 0);
! 139: size -= (fp[size-1] == 0);
! 140: SIZ(fn) = size;
1.1 maekawa 141:
1.1.1.2 ! ohara 142: TRACE (printf ("done special, size=%ld\n", size);
! 143: mpn_trace ("fp ", fp, size));
1.1 maekawa 144:
1.1.1.2 ! ohara 145: TMP_FREE (marker);
1.1 maekawa 146: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to fib_ui.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz