From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- gmp-6.3.0/mpn/generic/strongfibo.c | 219 +++++++++++++++++++++++++++++++++++++ 1 file changed, 219 insertions(+) create mode 100644 gmp-6.3.0/mpn/generic/strongfibo.c (limited to 'gmp-6.3.0/mpn/generic/strongfibo.c') diff --git a/gmp-6.3.0/mpn/generic/strongfibo.c b/gmp-6.3.0/mpn/generic/strongfibo.c new file mode 100644 index 0000000..7e8d612 --- /dev/null +++ b/gmp-6.3.0/mpn/generic/strongfibo.c @@ -0,0 +1,219 @@ +/* mpn_fib2m -- calculate Fibonacci numbers, modulo m. + +Contributed to the GNU project by Marco Bodrato. + + THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST + CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN + FUTURE GNU MP RELEASES. + +Copyright 2001, 2002, 2005, 2009, 2018, 2022 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include +#include "gmp-impl.h" + + +#if ! HAVE_NATIVE_mpn_rsblsh1_n && ! HAVE_NATIVE_mpn_sublsh1_n +/* Stores |{ap,n}-{bp,n}| in {rp,n}, + returns the sign of {ap,n}-{bp,n}. */ +static int +abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n) +{ + mp_limb_t x, y; + while (--n >= 0) + { + x = ap[n]; + y = bp[n]; + if (x != y) + { + ++n; + if (x > y) + { + ASSERT_NOCARRY (mpn_sub_n (rp, ap, bp, n)); + return 1; + } + else + { + ASSERT_NOCARRY (mpn_sub_n (rp, bp, ap, n)); + return -1; + } + } + rp[n] = 0; + } + return 0; +} +#endif + +/* Computes at most count terms of the sequence needed by the + Lucas-Lehmer-Riesel test, indexing backward: + L_i = L_{i+1}^2 - 2 + + The sequence is computed modulo M = {mp, mn}. + The starting point is given in L_{count+1} = {lp, mn}. + The scratch pointed by sp, needs a space of at least 3 * mn + 1 limbs. + + Returns the index i>0 if L_i = 0 (mod M) is found within the + computed count terms of the sequence. Otherwise it returns zero. + + Note: (+/-2)^2-2=2, (+/-1)^2-2=-1, 0^2-2=-2 + */ + +static mp_bitcnt_t +mpn_llriter (mp_ptr lp, mp_srcptr mp, mp_size_t mn, mp_bitcnt_t count, mp_ptr sp) +{ + do + { + mpn_sqr (sp, lp, mn); + mpn_tdiv_qr (sp + 2 * mn, lp, 0, sp, 2 * mn, mp, mn); + if (lp[0] < 5) + { + /* If L^2 % M < 5, |L^2 % M - 2| <= 2 */ + if (mn == 1 || mpn_zero_p (lp + 1, mn - 1)) + return (lp[0] == 2) ? count : 0; + else + MPN_DECR_U (lp, mn, 2); + } + else + lp[0] -= 2; + } while (--count != 0); + return 0; +} + +/* Store the Lucas' number L[n] at lp (maybe), computed modulo m. lp + and scratch should have room for mn*2+1 limbs. + + Returns the size of L[n] normally. + + If F[n] is zero modulo m, or L[n] is, returns 0 and lp is + undefined. +*/ + +static mp_size_t +mpn_lucm (mp_ptr lp, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn, mp_ptr scratch) +{ + int neg; + mp_limb_t cy; + + ASSERT (! MPN_OVERLAP_P (lp, MAX(2*mn+1,5), scratch, MAX(2*mn+1,5))); + ASSERT (nn > 0); + + neg = mpn_fib2m (lp, scratch, np, nn, mp, mn); + + /* F[n] = +/-{lp, mn}, F[n-1] = +/-{scratch, mn} */ + if (mpn_zero_p (lp, mn)) + return 0; + + if (neg) /* One sign is opposite, use sub instead of add. */ + { +#if HAVE_NATIVE_mpn_rsblsh1_n || HAVE_NATIVE_mpn_sublsh1_n +#if HAVE_NATIVE_mpn_rsblsh1_n + cy = mpn_rsblsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */ +#else + cy = mpn_sublsh1_n (lp, lp, scratch, mn); /* L[n] = -/+(F[n]-(-2F[n-1])) */ + if (cy != 0) + cy = mpn_add_n (lp, lp, mp, mn) - cy; +#endif + if (cy > 1) + cy += mpn_add_n (lp, lp, mp, mn); +#else + cy = mpn_lshift (scratch, scratch, mn, 1); /* 2F[n-1] */ + if (UNLIKELY (cy)) + cy -= mpn_sub_n (lp, scratch, lp, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */ + else + abs_sub_n (lp, lp, scratch, mn); +#endif + ASSERT (cy <= 1); + } + else + { +#if HAVE_NATIVE_mpn_addlsh1_n + cy = mpn_addlsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]+F[n])) */ +#else + cy = mpn_lshift (scratch, scratch, mn, 1); + cy+= mpn_add_n (lp, lp, scratch, mn); +#endif + ASSERT (cy <= 2); + } + while (cy || mpn_cmp (lp, mp, mn) >= 0) + cy -= mpn_sub_n (lp, lp, mp, mn); + MPN_NORMALIZE (lp, mn); + return mn; +} + +int +mpn_strongfibo (mp_srcptr mp, mp_size_t mn, mp_ptr scratch) +{ + mp_ptr lp, sp; + mp_size_t en; + mp_bitcnt_t b0; + TMP_DECL; + +#if GMP_NUMB_BITS % 4 == 0 + b0 = mpn_scan0 (mp, 0); +#else + { + mpz_t m = MPZ_ROINIT_N(mp, mn); + b0 = mpz_scan0 (m, 0); + } + if (UNLIKELY (b0 == mn * GMP_NUMB_BITS)) + { + en = 1; + scratch [0] = 1; + } + else +#endif + { + int cnt = b0 % GMP_NUMB_BITS; + en = b0 / GMP_NUMB_BITS; + if (LIKELY (cnt != 0)) + mpn_rshift (scratch, mp + en, mn - en, cnt); + else + MPN_COPY (scratch, mp + en, mn - en); + en = mn - en; + scratch [0] |= 1; + en -= scratch [en - 1] == 0; + } + TMP_MARK; + + lp = TMP_ALLOC_LIMBS (4 * mn + 6); + sp = lp + 2 * mn + 3; + en = mpn_lucm (sp, scratch, en, mp, mn, lp); + if (en != 0 && LIKELY (--b0 != 0)) + { + mpn_sqr (lp, sp, en); + lp [0] |= 2; /* V^2 + 2 */ + if (LIKELY (2 * en >= mn)) + mpn_tdiv_qr (sp, lp, 0, lp, 2 * en, mp, mn); + else + MPN_ZERO (lp + 2 * en, mn - 2 * en); + if (! mpn_zero_p (lp, mn) && LIKELY (--b0 != 0)) + b0 = mpn_llriter (lp, mp, mn, b0, lp + mn + 1); + } + TMP_FREE; + return (b0 != 0); +} -- cgit v1.2.3