From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001
From: Duncan Wilkie <antigravityd@gmail.com>
Date: Sat, 18 Nov 2023 06:11:09 -0600
Subject: Initial commit.

---
 gmp-6.3.0/mpz/fac_ui.c | 122 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 122 insertions(+)
 create mode 100644 gmp-6.3.0/mpz/fac_ui.c

(limited to 'gmp-6.3.0/mpz/fac_ui.c')

diff --git a/gmp-6.3.0/mpz/fac_ui.c b/gmp-6.3.0/mpz/fac_ui.c
new file mode 100644
index 0000000..a7998b3
--- /dev/null
+++ b/gmp-6.3.0/mpz/fac_ui.c
@@ -0,0 +1,122 @@
+/* mpz_fac_ui(RESULT, N) -- Set RESULT to N!.
+
+Contributed to the GNU project by Marco Bodrato.
+
+Copyright 1991, 1993-1995, 2000-2003, 2011, 2012, 2015, 2021 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 "gmp-impl.h"
+
+#define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I)		\
+  do {								\
+    if ((PR) > (MAX_PR)) {					\
+      (VEC)[(I)++] = (PR);					\
+      (PR) = (P);						\
+    } else							\
+      (PR) *= (P);						\
+  } while (0)
+
+#if TUNE_PROGRAM_BUILD
+#define FACTORS_PER_LIMB (GMP_NUMB_BITS / (LOG2C(FAC_DSC_THRESHOLD_LIMIT-1)+1))
+#else
+#define FACTORS_PER_LIMB (GMP_NUMB_BITS / (LOG2C(FAC_ODD_THRESHOLD)+1))
+#endif
+
+/* Computes n!, the factorial of n.
+   WARNING: it assumes that n fits in a limb!
+ */
+void
+mpz_fac_ui (mpz_ptr x, unsigned long n)
+{
+  static const mp_limb_t table[] = { ONE_LIMB_FACTORIAL_TABLE };
+
+  ASSERT (n <= GMP_NUMB_MAX);
+
+  if (n < numberof (table))
+    {
+      MPZ_NEWALLOC (x, 1)[0] = table[n];
+      SIZ (x) = 1;
+    }
+  else if (BELOW_THRESHOLD (n, FAC_ODD_THRESHOLD))
+    {
+      mp_limb_t prod, max_prod;
+      mp_size_t j;
+      mp_ptr    factors;
+      mp_limb_t fac, diff = n - numberof (table);
+      TMP_SDECL;
+
+      TMP_SMARK;
+      factors = TMP_SALLOC_LIMBS (2 + diff / FACTORS_PER_LIMB);
+
+      factors[0] = table[numberof (table)-1];
+      j = 1;
+      if ((diff & 1) == 0)
+	{
+	  prod = n;
+	  /* if (diff != 0) */
+	    fac = --n * numberof (table);
+	}
+      else
+	{
+	  prod = n * numberof (table);
+	  fac = prod + --diff;
+	}
+
+#if TUNE_PROGRAM_BUILD
+      max_prod = GMP_NUMB_MAX / (FAC_DSC_THRESHOLD_LIMIT * FAC_DSC_THRESHOLD_LIMIT);
+#else
+      max_prod = GMP_NUMB_MAX /
+	(((FAC_ODD_THRESHOLD + numberof (table) + 1) / 2) *
+	 ((FAC_ODD_THRESHOLD + numberof (table)) / 2));
+#endif
+      for (;diff != 0; fac += (diff -= 2))
+	FACTOR_LIST_STORE (fac, prod, max_prod, factors, j);
+
+      factors[j++] = prod;
+      mpz_prodlimbs (x, factors, j);
+
+      TMP_SFREE;
+    }
+  else
+    {
+      mp_limb_t count;
+      mpz_oddfac_1 (x, n, 0);
+      if (n <= TABLE_LIMIT_2N_MINUS_POPC_2N)
+	count = __gmp_fac2cnt_table[n / 2 - 1];
+      else
+	{
+	  popc_limb (count, n);
+	  count = n - count;
+	}
+      mpz_mul_2exp (x, x, count);
+    }
+}
+
+#undef FACTORS_PER_LIMB
+#undef FACTOR_LIST_STORE
-- 
cgit v1.2.3