More accurate calculation of sieve array size

[Adrien Prost-Boucle]

Hi,

Looking into oddfac and primorial computation code,
I think there is a suboptimality in the computation of the prime sieve size.
Could anyone confirm?

Currently in promorial, the number of words allocated is basically this formula :
  n / GMP_NUMB_BITS * 1.5
But the sieve requires only this number of words :
  n / GMP_NUMB_BITS / 3

This is also what is done in nextprime.c, and in bin_uiui.c
The suboptimality also exists in oddfac.
Which leads to the proposed patch below.

Regards,
Adrien


--- mpz/oddfac_1.old	2022-08-10 15:14:04.511298108 +0200
+++ mpz/oddfac_1.c	2022-08-10 15:31:34.721354647 +0200
@@ -381,8 +381,8 @@
 	  TMP_MARK;
 
 	  flag--;
-	  size = n / GMP_NUMB_BITS + 4;
-	  ASSERT (primesieve_size (n - 1) <= size - (size / 2 + 1));
+	  size = (n / 3 / GMP_NUMB_BITS + 1) * 2;
+	  ASSERT (primesieve_size (n - 1) <= size / 2);
 	  /* 2-multiswing(n) < 2^(n-1)*sqrt(n/pi) < 2^(n+GMP_NUMB_BITS);
 	     one more can be overwritten by mul, another for the sieve */
 	  MPZ_TMP_INIT (mswing, size);
@@ -390,7 +390,7 @@
 	  ASSERT_CODE (SIZ (mswing) = 0);
 
 	  /* Put the sieve on the second half, it will be overwritten by the last mswing. */
-	  sieve = PTR (mswing) + size / 2 + 1;
+	  sieve = PTR (mswing) + size / 2;
 
 	  size = (gmp_primesieve (sieve, n - 1) + 1) / log_n_max (n) + 1;
 

--- mpz/primorial_ui.c.old	2022-08-10 14:44:04.437867901 +0200
+++ mpz/primorial_ui.c	2022-08-10 14:55:41.721238758 +0200
@@ -82,8 +82,7 @@
       mp_limb_t prod;
       TMP_DECL;
 
-      size = n / GMP_NUMB_BITS;
-      size = size + (size >> 1) + 1;
+      size = n / 3 / GMP_NUMB_BITS + 1;
       ASSERT (size >= primesieve_size (n));
       sieve = MPZ_NEWALLOC (x, size);
       size = (gmp_primesieve (sieve, n) + 1) / log_n_max (n) + 1;