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 } 168