1 /* mpn_fib2m -- calculate Fibonacci numbers, modulo m. 2 3 Contributed to the GNU project by Marco Bodrato. 4 5 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST 6 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN 7 FUTURE GNU MP RELEASES. 8 9 Copyright 2001, 2002, 2005, 2009, 2018 Free Software Foundation, Inc. 10 11 This file is part of the GNU MP Library. 12 13 The GNU MP Library is free software; you can redistribute it and/or modify 14 it under the terms of either: 15 16 * the GNU Lesser General Public License as published by the Free 17 Software Foundation; either version 3 of the License, or (at your 18 option) any later version. 19 20 or 21 22 * the GNU General Public License as published by the Free Software 23 Foundation; either version 2 of the License, or (at your option) any 24 later version. 25 26 or both in parallel, as here. 27 28 The GNU MP Library is distributed in the hope that it will be useful, but 29 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 30 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 31 for more details. 32 33 You should have received copies of the GNU General Public License and the 34 GNU Lesser General Public License along with the GNU MP Library. If not, 35 see https://www.gnu.org/licenses/. */ 36 37 #include <stdio.h> 38 #include "gmp-impl.h" 39 40 41 #if HAVE_NATIVE_mpn_rsblsh1_n || HAVE_NATIVE_mpn_sublsh1_n 42 #else 43 /* Stores |{ap,n}-{bp,n}| in {rp,n}, 44 returns the sign of {ap,n}-{bp,n}. */ 45 static int 46 abs_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 #endif 72 73 /* Computes at most count terms of the sequence needed by the 74 Lucas-Lehmer-Riesel test, indexing backward: 75 L_i = L_{i+1}^2 - 2 76 77 The sequence is computed modulo M = {mp, mn}. 78 The starting point is given in L_{count+1} = {lp, mn}. 79 The scratch pointed by sp, needs a space of at least 3 * mn + 1 limbs. 80 81 Returns the index i>0 if L_i = 0 (mod M) is found within the 82 computed count terms of the sequence. Otherwise it returns zero. 83 84 Note: (+/-2)^2-2=2, (+/-1)^2-2=-1, 0^2-2=-2 85 */ 86 87 static mp_bitcnt_t 88 mpn_llriter (mp_ptr lp, mp_srcptr mp, mp_size_t mn, mp_bitcnt_t count, mp_ptr sp) 89 { 90 do 91 { 92 mpn_sqr (sp, lp, mn); 93 mpn_tdiv_qr (sp + 2 * mn, lp, 0, sp, 2 * mn, mp, mn); 94 if (lp[0] < 5) 95 { 96 /* If L^2 % M < 5, |L^2 % M - 2| <= 2 */ 97 if (mn == 1 || mpn_zero_p (lp + 1, mn - 1)) 98 return (lp[0] == 2) ? count : 0; 99 else 100 MPN_DECR_U (lp, mn, 2); 101 } 102 else 103 lp[0] -= 2; 104 } while (--count != 0); 105 return 0; 106 } 107 108 /* Store the Lucas' number L[n] at lp (maybe), computed modulo m. lp 109 and scratch should have room for mn*2+1 limbs. 110 111 Returns the size of L[n] normally. 112 113 If F[n] is zero modulo m, or L[n] is, returns 0 and lp is 114 undefined. 115 */ 116 117 static mp_size_t 118 mpn_lucm (mp_ptr lp, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn, mp_ptr scratch) 119 { 120 int neg; 121 mp_limb_t cy; 122 123 ASSERT (! MPN_OVERLAP_P (lp, MAX(2*mn+1,5), scratch, MAX(2*mn+1,5))); 124 ASSERT (nn > 0); 125 126 neg = mpn_fib2m (lp, scratch, np, nn, mp, mn); 127 128 /* F[n] = +/-{lp, mn}, F[n-1] = +/-{scratch, mn} */ 129 if (mpn_zero_p (lp, mn)) 130 return 0; 131 132 if (neg) /* One sign is opposite, use sub instead of add. */ 133 { 134 #if HAVE_NATIVE_mpn_rsblsh1_n || HAVE_NATIVE_mpn_sublsh1_n 135 #if HAVE_NATIVE_mpn_rsblsh1_n 136 cy = mpn_rsblsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */ 137 #else 138 cy = mpn_sublsh1_n (lp, lp, scratch, mn); /* L[n] = -/+(F[n]-(-2F[n-1])) */ 139 if (cy != 0) 140 cy = mpn_add_n (lp, lp, mp, mn) - cy; 141 #endif 142 if (cy > 1) 143 cy += mpn_add_n (lp, lp, mp, mn); 144 #else 145 cy = mpn_lshift (scratch, scratch, mn, 1); /* 2F[n-1] */ 146 if (UNLIKELY (cy)) 147 cy -= mpn_sub_n (lp, scratch, lp, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */ 148 else 149 abs_sub_n (lp, lp, scratch, mn); 150 #endif 151 ASSERT (cy <= 1); 152 } 153 else 154 { 155 #if HAVE_NATIVE_mpn_addlsh1_n 156 cy = mpn_addlsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]+F[n])) */ 157 #else 158 cy = mpn_lshift (scratch, scratch, mn, 1); 159 cy+= mpn_add_n (lp, lp, scratch, mn); 160 #endif 161 ASSERT (cy <= 2); 162 } 163 while (cy || mpn_cmp (lp, mp, mn) >= 0) 164 cy -= mpn_sub_n (lp, lp, mp, mn); 165 MPN_NORMALIZE (lp, mn); 166 return mn; 167 } 168 169 int 170 mpn_strongfibo (mp_srcptr mp, mp_size_t mn, mp_ptr scratch) 171 { 172 mp_ptr lp, sp; 173 mp_size_t en; 174 mp_bitcnt_t b0; 175 TMP_DECL; 176 177 #if GMP_NUMB_BITS % 4 == 0 178 b0 = mpn_scan0 (mp, 0); 179 #else 180 { 181 mpz_t m = MPZ_ROINIT_N(mp, mn); 182 b0 = mpz_scan0 (m, 0); 183 } 184 if (UNLIKELY (b0 == mn * GMP_NUMB_BITS)) 185 { 186 en = 1; 187 scratch [0] = 1; 188 } 189 else 190 #endif 191 { 192 int cnt = b0 % GMP_NUMB_BITS; 193 en = b0 / GMP_NUMB_BITS; 194 if (LIKELY (cnt != 0)) 195 mpn_rshift (scratch, mp + en, mn - en, cnt); 196 else 197 MPN_COPY (scratch, mp + en, mn - en); 198 en = mn - en; 199 scratch [0] |= 1; 200 en -= scratch [en - 1] == 0; 201 } 202 TMP_MARK; 203 204 lp = TMP_ALLOC_LIMBS (4 * mn + 6); 205 sp = lp + 2 * mn + 3; 206 en = mpn_lucm (sp, scratch, en, mp, mn, lp); 207 if (en != 0 && LIKELY (--b0 != 0)) 208 { 209 mpn_sqr (lp, sp, en); 210 lp [0] |= 2; /* V^2 + 2 */ 211 if (LIKELY (2 * en >= mn)) 212 mpn_tdiv_qr (sp, lp, 0, lp, 2 * en, mp, mn); 213 else 214 MPN_ZERO (lp + 2 * en, mn - 2 * en); 215 if (! mpn_zero_p (lp, mn) && LIKELY (--b0 != 0)) 216 b0 = mpn_llriter (lp, mp, mn, b0, lp + mn + 1); 217 } 218 TMP_FREE; 219 return (b0 != 0); 220 } 221
Indexes created Sat Feb 28 05:31:39 UTC 2026