1/* mpn_fib2m -- calculate Fibonacci numbers, modulo m. 2 3Contributed to the GNU project by Marco Bodrato, based on the previous 4fib2_ui.c file. 5 6 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST 7 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN 8 FUTURE GNU MP RELEASES. 9 10Copyright 2001, 2002, 2005, 2009, 2018 Free Software Foundation, Inc. 11 12This file is part of the GNU MP Library. 13 14The GNU MP Library is free software; you can redistribute it and/or modify 15it under the terms of either: 16 17 * the GNU Lesser General Public License as published by the Free 18 Software Foundation; either version 3 of the License, or (at your 19 option) any later version. 20 21or 22 23 * the GNU General Public License as published by the Free Software 24 Foundation; either version 2 of the License, or (at your option) any 25 later version. 26 27or both in parallel, as here. 28 29The GNU MP Library is distributed in the hope that it will be useful, but 30WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 31or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 32for more details. 33 34You should have received copies of the GNU General Public License and the 35GNU Lesser General Public License along with the GNU MP Library. If not, 36see https://www.gnu.org/licenses/. */ 37 38#include <stdio.h> 39#include "gmp-impl.h" 40#include "longlong.h" 41 42 43/* Stores |{ap,n}-{bp,n}| in {rp,n}, 44 returns the sign of {ap,n}-{bp,n}. */ 45static int 46abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) 47{ 48 mp_limb_t x, y; 49 while (--n >= 0) 50 { 51 x = ap[n]; 52 y = bp[n]; 53 if (x != y) 54 { 55 ++n; 56 if (x > y) 57 { 58 ASSERT_NOCARRY (mpn_sub_n (rp, ap, bp, n)); 59 return 1; 60 } 61 else 62 { 63 ASSERT_NOCARRY (mpn_sub_n (rp, bp, ap, n)); 64 return -1; 65 } 66 } 67 rp[n] = 0; 68 } 69 return 0; 70} 71 72/* Store F[n] at fp and F[n-1] at f1p. Both are computed modulo m. 73 fp and f1p should have room for mn*2+1 limbs. 74 75 The sign of one or both the values may be flipped (n-F, instead of F), 76 the return value is 0 (zero) if the signs are coherent (both positive 77 or both negative) and 1 (one) otherwise. 78 79 Notes: 80 81 In F[2k+1] with k even, +2 is applied to 4*F[k]^2 just by ORing into the 82 low limb. 83 84 In F[2k+1] with k odd, -2 is applied to F[k-1]^2 just by ORing into the 85 low limb. 86 87 TODO: Should {tp, 2 * mn} be passed as a scratch pointer? 88 Should the call to mpn_fib2_ui() obtain (up to) 2*mn limbs? 89*/ 90 91int 92mpn_fib2m (mp_ptr fp, mp_ptr f1p, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn) 93{ 94 unsigned long nfirst; 95 mp_limb_t nh; 96 mp_bitcnt_t nbi; 97 mp_size_t sn, fn; 98 int fcnt, ncnt; 99 100 ASSERT (! MPN_OVERLAP_P (fp, MAX(2*mn+1,5), f1p, MAX(2*mn+1,5))); 101 ASSERT (nn > 0 && np[nn - 1] != 0); 102 103 /* Estimate the maximal n such that fibonacci(n) fits in mn limbs. */ 104#if GMP_NUMB_BITS % 16 == 0 105 if (UNLIKELY (ULONG_MAX / (23 * (GMP_NUMB_BITS / 16)) <= mn)) 106 nfirst = ULONG_MAX; 107 else 108 nfirst = mn * (23 * (GMP_NUMB_BITS / 16)); 109#else 110 { 111 mp_bitcnt_t mbi; 112 mbi = (mp_bitcnt_t) mn * GMP_NUMB_BITS; 113 114 if (UNLIKELY (ULONG_MAX / 23 < mbi)) 115 { 116 if (UNLIKELY (ULONG_MAX / 23 * 16 <= mbi)) 117 nfirst = ULONG_MAX; 118 else 119 nfirst = mbi / 16 * 23; 120 } 121 else 122 nfirst = mbi * 23 / 16; 123 } 124#endif 125 126 sn = nn - 1; 127 nh = np[sn]; 128 count_leading_zeros (ncnt, nh); 129 count_leading_zeros (fcnt, nfirst); 130 131 if (fcnt >= ncnt) 132 { 133 ncnt = fcnt - ncnt; 134 nh >>= ncnt; 135 } 136 else if (sn > 0) 137 { 138 ncnt -= fcnt; 139 nh <<= ncnt; 140 ncnt = GMP_NUMB_BITS - ncnt; 141 --sn; 142 nh |= np[sn] >> ncnt; 143 } 144 else 145 ncnt = 0; 146 147 nbi = sn * GMP_NUMB_BITS + ncnt; 148 if (nh > nfirst) 149 { 150 nh >>= 1; 151 ++nbi; 152 } 153 154 ASSERT (nh <= nfirst); 155 /* Take a starting pair from mpn_fib2_ui. */ 156 fn = mpn_fib2_ui (fp, f1p, nh); 157 MPN_ZERO (fp + fn, mn - fn); 158 MPN_ZERO (f1p + fn, mn - fn); 159 160 if (nbi == 0) 161 { 162 if (fn == mn) 163 { 164 mp_limb_t qp[2]; 165 mpn_tdiv_qr (qp, fp, 0, fp, fn, mp, mn); 166 mpn_tdiv_qr (qp, f1p, 0, f1p, fn, mp, mn); 167 } 168 169 return 0; 170 } 171 else 172 { 173 mp_ptr tp; 174 unsigned pb = nh & 1; 175 int neg; 176 TMP_DECL; 177 178 TMP_MARK; 179 180 tp = TMP_ALLOC_LIMBS (2 * mn + (mn < 2)); 181 182 do 183 { 184 mp_ptr rp; 185 /* Here fp==F[k] and f1p==F[k-1], with k being the bits of n from 186 nbi upwards. 187 188 Based on the next bit of n, we'll double to the pair 189 fp==F[2k],f1p==F[2k-1] or fp==F[2k+1],f1p==F[2k], according as 190 that bit is 0 or 1 respectively. */ 191 192 mpn_sqr (tp, fp, mn); 193 mpn_sqr (fp, f1p, mn); 194 195 /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2. */ 196 f1p[2 * mn] = mpn_add_n (f1p, tp, fp, 2 * mn); 197 198 /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k. 199 pb is the low bit of our implied k. */ 200 201 /* fp is F[k-1]^2 == 0 or 1 mod 4, like all squares. */ 202 ASSERT ((fp[0] & 2) == 0); 203 ASSERT (pb == (pb & 1)); 204 ASSERT ((fp[0] + (pb ? 2 : 0)) == (fp[0] | (pb << 1))); 205 fp[0] |= pb << 1; /* possible -2 */ 206#if HAVE_NATIVE_mpn_rsblsh2_n 207 fp[2 * mn] = 1 + mpn_rsblsh2_n (fp, fp, tp, 2 * mn); 208 MPN_INCR_U(fp, 2 * mn + 1, (1 ^ pb) << 1); /* possible +2 */ 209 fp[2 * mn] = (fp[2 * mn] - 1) & GMP_NUMB_MAX; 210#else 211 { 212 mp_limb_t c; 213 214 c = mpn_lshift (tp, tp, 2 * mn, 2); 215 tp[0] |= (1 ^ pb) << 1; /* possible +2 */ 216 c -= mpn_sub_n (fp, tp, fp, 2 * mn); 217 fp[2 * mn] = c & GMP_NUMB_MAX; 218 } 219#endif 220 neg = fp[2 * mn] == GMP_NUMB_MAX; 221 222 /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2 */ 223 /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k */ 224 225 /* Calculate F[2k] = F[2k+1] - F[2k-1], replacing the unwanted one of 226 F[2k+1] and F[2k-1]. */ 227 --nbi; 228 pb = (np [nbi / GMP_NUMB_BITS] >> (nbi % GMP_NUMB_BITS)) & 1; 229 rp = pb ? f1p : fp; 230 if (neg) 231 { 232 /* Calculate -(F[2k+1] - F[2k-1]) */ 233 rp[2 * mn] = f1p[2 * mn] + 1 - mpn_sub_n (rp, f1p, fp, 2 * mn); 234 neg = ! pb; 235 if (pb) /* fp not overwritten, negate it. */ 236 fp [2 * mn] = 1 ^ mpn_neg (fp, fp, 2 * mn); 237 } 238 else 239 { 240 neg = abs_sub_n (rp, fp, f1p, 2 * mn + 1) < 0; 241 } 242 243 mpn_tdiv_qr (tp, fp, 0, fp, 2 * mn + 1, mp, mn); 244 mpn_tdiv_qr (tp, f1p, 0, f1p, 2 * mn + 1, mp, mn); 245 } 246 while (nbi != 0); 247 248 TMP_FREE; 249 250 return neg; 251 } 252} 253