Initial commit.

1 files changed, 219 insertions, 0 deletions

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 <stdio.h>
+#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);
+}